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Abstract 



In this paper we prove that parabolic Julia sets of rational functions are locally 
computable in polynomial time. 

1 Introduction 

In the present paper we consider the complexity of generating precise images of Julia sets 
with parabolic orbits. It has been independently proved in [Brv04j and |Ret04| that hyperbolic 
Julia sets can be computed in polynomial time. Neither of the two algorithms can be applied 
in the parabolic case. In fact, both algorithms often slow down significantly as the underlying 
polynomial approaches one with a parabolic point. A naive generalization of these algorithms 
would yield exponential time algorithms in the parabolic case, which are useless when one 
is trying to produce meaningful pictures of the Julia set in question. 

The same problem has been highlighted in the comments on computer graphics by John 
Milnor in (MilOOj . Appendix H. The example considered there is for the polynomial p(z) = 
z + z A . It has a parabolic fixed point at z — 0. Consider a point z = e ~ 1/1000. Suppose 
we are trying to determine whether e is in the Julia set or not by iterating it, and observing 
whether its orbit escapes to oo, or converges to 0. In fact, such a z would always escape to 
oo, but it is not hard to see that this process would take l/3e 3 ~ 300,000,000 iterations 
for z to escape the ball of radius 2 around 0. Thus, we would need to follow the orbit 
for « 300, 000, 000 iterations before concluding that it converges to oo. If we zoom-in a 
little and set z = e ~ 1/100,000, we would need ~ 3 • 10 14 iterations to trace z, which is 
computationally impractical. 

Due to the effect highlighted above, most computer programs plotting Julia sets include 
all the points that diverge slowly from the parabolic orbit in the Julia set. 

The algorithm we present here is not uniform, i.e. it requires a special program for each 
specific parabolic Julia set. The running time of the algorithm is C r n c , where the constant 
C r depends on the rational function r but not on n, and c is some constant. The algorithm 
can be made uniform in the r, provided some basic combinatorial information about the 
parabolic points. I.e. one algorithm can compute all the parabolic sets, if it is provided with 
some basic information about the rational function. The constant C r in the running time 
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can still vary strongly for different functions r. For example, it is reasonable to expect that 
J n for ri(z) = z + z 2 would take less time to compute than J r2 for r 2 (z) = e 2m l 17 z + z 23 . 
We prove the following: 

Theorem 1 There is an algorithm A that given 

• a rational function r(z) such that every critical orbit of r converges either to an at- 
tracting or a parabolic orbit; and 

• some basic combinatorial information about the parabolic orbits of r; 

produces an image of the Julia set J r . A takes time C r n c to decide one pixel in J r with 
precision 2~ n . Here c is some small constant and C r depends on r but not on n. 

After this work was completed, John Milnor has informed us that he has used an algo- 
rithm similar to ours to produce pictures of Julia sets with parabolic points. In particular, 
some of the pictures in |MilO Ct were created this way. 

The rest of the paper is organized as follows. In section |21 we give the necessary prelimi- 
naries on the complexity theory over the reals. In section El we outline the general strategy 
for computing parabolic Julia sets fast. Sections IH El and El provide the main tool for the 
algorithm - computing a "long" iteration near a parabolic point. Finally, in sectional we 
present and analyze the algorithm. 

Acknowledgment. The author wishes to thank Ilia Binder and Michael Yampolsky for 
their insights and encouragement during the preparation of this paper. 

2 Complexity over R — preliminaries 

In this section we provide some preliminaries on the notion of complexity for sets and func- 
tions over R n , in particular R 2 . More details can be found in jB W99j . jBrv05j and [WeiOOj. 

2.1 Complexity of Sets in R 2 

Intuitively, we say the computational complexity of a set S is t(n) if it takes time t(n) to 
decide whether to draw a pixel of size 2~ n in the picture of S. To make this notion precise, 
we have to decide what are our expectations from a picture of S. First of all, we expect a 
good picture of S to cover the whole set S. On the other hand, we expect every point of 
the picture to be close to some point of S, otherwise the picture would have no descriptive 
power about S. Mathematically, we write these requirements as follows: 

Definition 2 A set T is said to be a 2~ n -picture of a bounded set S if 

(i) S C T, and (ii) T C B(S, 2~ n ) = {x e R 2 : \x - s\ < 2~ n for some s E S}. 
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Definition El is also equivalent to approximating S by T in the Hausdorff metric, given 

by 

d H (S,T) := inf {r : S C B (T,r) and Tcfi^r)}. 

Suppose we are trying to generate a picture of a set S 1 using a union of round pixels 
of radius 2~ n with centers at all the points of the form ^J, with i and j integers. In 
order to draw the picture, we have to decide for each pair whether to draw the pixel 
centered at (j^, -^j or not. We want to draw the pixel if it intersects S and to omit it if some 
neighborhood of the pixel does not intersect S. Formally, we want to compute a function 

f 1, 5((V2^/2"),2-")nSV0 
f s (n,z/2 n ,j/2 n ) = l 0, S((</2» j/2») ) 2.2-»)n5 = (1) 

| or 1, in all other cases 




Figure 1: Sample values of /. The radius of the inner circle is 2 n 2 . 

Lemma 3 The picture drawn according to fs{n, •) is a 2~( n ~ 2 ^ -picture of S. 

Here • stands for the different values of the parameters (i/2 n , j/2 n ). The lemma illustrates 
the tight connection between the complexity of "drawing" the set S and the complexity of 
computing /. We reflect this connection by defining the time complexity of S as follows. 

Definition 4 A bounded set S is said to be computable in time t(n) if there is a function 
f(n, •) satisfying (QJj which runs in time t(n). We say that S is poly-time computable if there 
is a polynomial p, such that S is computable in time p(n). 

To see why this is the "right" definition, suppose we are trying to draw a set S on a 
computer screen which has a 1000 x 1000 pixel resolution. A 2 -n -zoomed in picture of S has 
0(2 2n ) pixels of size 2~ n , and thus would take time 0(t(n) ■ 2~ 2n ) to compute. This quantity 
is exponential in n, even if t(n) is bounded by a polynomial. But we are drawing S on a 
finite-resolution display, and we will only need to draw 1000 • 1000 = 10 6 pixels. Hence the 
running; time would be O(10 6 • t{n)) = 0(t(n)). This running time is polynomial in n if and 
only if t(n) is polynomial. Hence t(n) reflects the 'true' cost of zooming in. 
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2.2 Computing Julia Sets 



There are uncountably many rational functions, but only countably many Turing Machines. 
Thus, we cannot expect to have a Turing Machine computing the Julia set J r for each rational 
r(z). Instead, we assume that the coefficients of r are given to the machine, and it is trying 
to produce a picture of J r . The machine can access the coefficients with an arbitrarily high 
(finite precision). It is charged m time units for querying a coefficient with precision 2~ m . 
Hence if a machine computes J r with precision 2~ n in time polynomial p(n), it will query 
the coefficients with precision at most 2~ p(jl \ 

Another issue is whether the computation of a machine is uniform or non-uniform. A 
machine for computing J r is non-uniform, if it is designed specifically for this r. A machine 
is uniform on a set S of rational functions, if it produces J r for all r G S. One can view a 
non- uniform machine as a uniform machine on the set S = {r}. One of the properties of the 
computation model is that if J r is uniformly computable on S, then the function J : r i— ► J r 
is continuous in the Hausdorff metric. In the case of a non-uniform computation, S is a 
singleton, and thus we don't get any information from this statement. 

We first give a non- uniform algorithm for computing J r . Then in section we argue that 
it can be made uniform for some large classes of parabolic Julia sets. The function J : r i— ► J r 
is not continuous over all parabolic sets, and thus it cannot be uniformly computable on all 
parabolic functions r. See section I77S1 for more details. 

3 The Strategy 

First we recall the strategy in the hyperbolic case, which is much easier to deal with. Suppose 
that r is a hyperbolic rational function. Let J r denote its Julia set. Then r is strictly 
expanding by some constant c > 1 in the hyperbolic metric around J r , and thus the escape 
rate of a point z £ J r near J r is exponential. In other words, if d(z, J r ) > 2~ n , then after 
0(n) steps the orbit of z will be at 0(1) distance from J r . This gives a natural poly-time 
algorithm for computing J r : iterate z until it is possible to estimate the distance from 
r k (z) to J r using some coarse initial approximation to J r . If such a k = 0(n) exists, use 
d(r k (z), J r ) and |(r fc )'(;z)| to estimate d(z, J r ). If no such k exists, we can be sure that initially 
d(z, J r ) < 2~ n . 

We would like to employ a similar strategy here, in the parabolic case. The problem is 
that even though r is still expanding in the hyperbolic metric near J r , the expansion is now 
extremely slow near the parabolic point. For example, let r(z) = z 2 + l/4, with the parabolic 
point p = 1/2. The picture of J r is presented of figure 121 If we set z = 1/2 + 2"™, it will take 
0{2 n ) steps before z escapes the unit disk. 

We solve this problem by approximating a "long" iteration of z in the neighborhood of 
a parabolic point fast. In the previous example "long" would mean 0(2"). 

On figure |2l we present the different regions which will appear in the algorithm. We list 
them below. 

• J r is the Julia set we are trying to compute. 
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Figure 2: A schematic image of the components in the algorithm. 



• U is some small fixed region around J r . All the points in U are much closer to J r 
than to the postcritical set. U is bounded away from J r , except for a finite number of 
parabolic and preparabolic touching points. On figure El the touching points are 1/2 
(the parabolic point) and —1/2 (first-order preparabolic). 

• E is the region around the parabolic points in which the "long" iteration is applicable. 
We also include in E preimages of this neighborhood around the preparabolic touching 
points of U and J r . 

• A C UUE is a collection of small wedges around the repelling directions. These wedges 
contain the portions of J r in the neighborhood of the corresponding parabolic/preparabolic 
points. 

• O is a tiny neighborhood around the touching points of U and J r . If the orbit of z falls 
into O we can be sure that z is close to J r because all of O is so close to J r . 

Now the algorithm works exactly as the one in the hyperbolic case: 

1. Iterate the orbit of z\ 

2. let w = r k (z) be the current iterate; 
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3. if w ^ U U E, we can estimate its distance from J r in 0(1) time; 

4. if w G U — E, just make one step w <— r(w); 

5. if w G O, output "u> close to J r "; 

6. if iu G £TlA in the neighborhood of a preparabolic point, just make one step w <— r(w); 

7. if w e E — A, we can estimate its distance from J r in 0(1) time; 

8. if w G Ed A near some parabolic point, apply linearly many "long" iterations to escape 
this region and get to step |U 

Step 0] can only be executed linearly many times, since dll is bounded away from J r 
outside of E and the expansion in the hyperbolic metric is bounded from below by some 
c > 1 on U — E. Thus the entire computation takes at most a quadratic number of steps to 
complete (at most linearly many executions of step El between two executions of stepEJ). 

Of course, this is only a sketch, and we need more precise procedures taking into account 
the finite precision of the computation etc. (e.g. we cannot just check whether w is in A or 
not.) In the next sections we will develop the tools for performing the "long" iteration near 
the parabolic points, before formally presenting the algorithm. 



4 Controlling coefficient growth 

The primary goal of this section is to prove the following lemma. 

Lemma 5 Let r > 1 be some integer. Set f(z) = z + z r+l + z r+2 + z r+3 + . . .. Then there 
is an explicit a such that the coefficients of the n-th iteration f n of f , 

f n (z) = z + a r z r+1 + a r+1 z r+2 + a r+2 z r+3 + 

satisfy 

a k < {an) k /\ 

One can take a = 2r 3 . 



We begin with a very simple proof in the case r = 1. The general case is more involved. 

Proof: (in the case r = 1.) In this case f(z) = z + z 2 + z 3 + . . . = (within the region 
of convergence) . It is easy to verify that 



Hence the coefficient of z k+1 in f n is n k , and lemma El holds with a — 1. I 

For the rest of the section we fix some r > 2, for which we are proving lemma El We will 
prove the lemma by induction on n. It is obviously true for n = 1 with any a > 1. Denote 
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p = (an) l l r . We will find a constant a that works later in the proof, a may depend on r but 
not on n. 

If g(z) and h(z) are two power series with positive real coefficients, we say that g is 
dominated by h, and write g -C h if all the coefficients of g are smaller or equal to the 
corresponding coefficient of h. 

We assume by the induction hypothesis that 

f n (z)^z + p r z r+1 +p r+1 z r+2 + .... (2) 

Denote g(z) = z + rz r+1 + r2 2r+1 + rz 3r+1 + . . .. We claim the following. 

Lemma 6 Let m > r be a given integer number. Then 

f(z) m - z m < (1 + z + z 2 + . . . + z r ~ l ) ■ [g{z) m - z m ). 

Proof: We show that the coefficient of z k on the left hand side is smaller or equal to the 
coefficient on the right hand side. Note that all the coefficients on the left for k < m + r are 
0, and so we can assume that k > m + r. Write k — m = r ■ I + q with < q < r, I > 1. We 
claim that the coefficient of z k ~ q = z r ' l+m in g(z) m is bigger or equal to the coefficient of z k 
in f(z) m . 

To see this we create a one-to-one mapping from all the terms of degree k in the expansion 
of f{z) m to the corresponding terms of degree k — q in the expansion of g(z) k where we write 

n (~\ - 7 1 i 7 r+l , r+1 , , r+1 , 2r+l , 2r+l , , 2r+l , 
fi^J - 2 (0) + Z (0) + Z (l) + • • • + Z (r-l) + Z (0) + 2 (1) + • • • + 2 (r-l) + 

Here Z( ), . . . , -Z( r -i) refer to different copies of the same z (we separate the different copies 
to specify the one-to-one mapping). 

Suppose we are given a term z ai z a2 . . . z am = z k in the expansion of f(z) m . We write 
dj = 1 + r • hi + Cj, < Cj < r. Then we know that either a« — q = 1 and c« = 0, or 
o>i — Ci > r + 1. We associate the term 

a\— ci _02— C2 Q m-i— c m _i a m +ci+C2+...+c fc _!— q /r>\ 

^(C2) ■•■ Z ( Cm -i) Z (0) ■ 

By the construction a m + C\ + c 2 + . . . + c m _i — q = a m + ai + a 2 + . . . + a m _i — (m — 1) — q = 
k — m — q + 1 = 1 (mod r). k > m + r, so we will never need the term z m from g(z) m . 
It is not hard to see that the correspondence is one-to-one, since the information in (jSJ) is 
sufficient to recover the values of a±, . . . , a m . 

By considering the term z q ■ g(z) m , we complete the proof. I 

We are now ready to make the induction step in lemma El By the induction hypothesis 
and lemma EJ we have 

f n+ \z) < f(z) + p r z r+1 + p r {f(z) r+1 - z r+l ) + p r+1 z r+2 + p r+ \f(z) r+2 - z r+2 ) + . . . < 

f(z) + p r z r+1 + p r+1 z r+2 + p r+2 z r+3 + ...+ 
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+ (1 + z + z 2 + . . . + z r ~ x ) ■ [p r (g(z) r+1 - z r+1 ) + p r+1 (g(z) r+2 - z r+I ) + . . . 

Our goal is to bound the coefficient of z k , k > r + 1, in f n+1 (z) = f n (f(z)). The 
contribution from f(z) is always 1. We consider the contribution from g(z) m — z m , r + 1 < 
m < k. Write k = m + r- l + q, < q < r. Then we must have z q in the product, and the 
coefficient is the coefficient of z m+r ' 1 in p m ~ 1 g(z) m , which is bounded by p m_1 ( m+ / -1 ) r '- The 
contribution is nonzero only if / > 0. Thus, the coefficient is bounded by 

^ ^ k-q-r-l j 

1 + + J2 £ P k - rl - q -V fk-rl-q + l-l\ < 

9=0 2=1 V ' / 

^ 1 +^ 1 |W(*7 1 ) =^§>^C 7 1 ) ^ 

To prove the lemma, we need the condition 

(p+^l) <W« + l)) (M/r (4) 
Recall that p = (an) 1//r , hence we can rewrite (@J) as 

H 1/r + (^w)^ a '" + 1 )- 



We have 



r 2 \ / r 2 \ r r 3 / 2r 3 ' 
(an) 1 ^ + t — -7 — TT-r ) = (an) ■ 1 H ) < (an) • < (an) • 1 H 



an 



the last inequality holds whenever — < In 2. Finally, if we take a > 2r 3 , then 



(an) ■ f 1 + — | < (an) • ( 1 + = a(n + 1) 
\ an / \ nj 



as required. 



5 Computing the n-th iteration of / 

Suppose that we are given a function / presented as a power series, finite or infinite, f(z) = 
z + a 2 z 2 + a 3 z 3 + . . ., d > 2. Denote for the n-th iteration f n of / 

f n (z) = z + a% l) z 2 + a ( £ l) z 3 + ... (5) 



8 



The goal of this section is to show how to compute the values of ajj. with a given precision 
2~ l fast - in time polynomial in k, I, and logn. We will need this in order to interpolate 
"long" iterations of / around the parabolic point (0 in this case). 

First, we show that if / has a non-negative radius of convergence R, then we can assume 
that \cii\ < 1 for all i, with a fairly small overhead. We know that Yli a i(R/^Y converges. 
Hence there is a bound B > 1 such that a^R/2) 1 < B for all i. In other words, |ai| < 
B ■ (2/Rf < (2B/R) l /B. In case / is a rational function, it is easy to approximate the 
number 2B/R, or some power of two C = 2 C , such that |aj| < C ,_1 for all i. Conjugate / by 
the map z i-> Cz to obtain g(z) = Cf(z/C). Then f(z) = g(Cz)/C, and g n (z) = Cf n (z/C). 
The Taylor expansion of g(z) is 

g(z) = z + ^z 2 + ^z 3 + ... 

We see that all the coefficients of g(z) do not exceed 1 in absolute value. The Taylor 
expansion of the n-th iteration of g is 

(n) (n) 

g n (z) = z + ^ + ^ + ... 

Thus, to compute aif 1 with precision 2~ l , we would need to compute the coefficient of z k in 
g n with precision 2~( l+c ' h \ This is a linear overhead, and if we can compute approximations 
for g in time poly (k, I, log n), we will also be able to do it for /. From now on, we assume 
that |a.j| < 1 for all i. 

We prove the following lemma. 

Lemma 7 Suppose f(z) = z + c^z 2 + a^z 3 + . . . is given by its power series. Then the 
coefficient as in |2P can be presented by a polynomial of the form 

4 n) = "o + a i n + a 2 n2 + ■■■ + at-in^ 1 , (6) 

and the values of a\ for j = 2, 3, . . . , k, i = 0, . . . , k — 1, can be computed with precision 2~ s 
in time polynomial in k and s. 

Proof: We prove the lemma by providing an iterative algorithm that computes the values 
of af. In the same time we prove that is indeed of the form as in equation On the 
j-th iteration we compute the values of a 3 , a{, . . . , 

Here is how to compute the a{ from the o4 _1 's. We know that 

riz) = f n -\f(z)) = fiz) + at 1] ■ f(zf + at 1] ■ f{z? + ■■■ 

In order to compute ai we need to find the coefficient of in each of the terms f(z), f(z) 2 , 
. . ., f(zY (it is always in higher terms). All we have to do is to compute all the coefficients 
of f(z), f(z) 2 , . . . , f{zy up to z^ with precision 2~ m in time polynomial in m and j. 
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We assume here that the coefficients a^-, a^, . . . , cij are given as oracles. In case that / is 
a rational function 

c — q{z) ■ z 

we know that c 7^ 0, by the parabolicity of / around 0, hence 

t/ \ 1 z + p(z)-z 2 1 ~ 2 /g(z) ■ z\ k 

c 1 - q{z) ■ z/c c ^ \ c J 

The first j coefficients of the expansion are now easily computed from this last formula. 

The computation is done using a simple "doubling" algorithm: first compute f(z) 2 , 
f(z) 22 , f(z) 23 , . . ., /(2;) 2L1 ° sjJ , and then compute the desired power of f(z) as a combination 
of these. At each multiplication we "chop" all the terms of degree j + 1 and above, hence 
the entire computation is polynomial. All the coefficients at all times are bounded by j- 7 in 
absolute value, hence the error is multiplied by at most 2ji = 2°^ log ^ at each step, and 
we will need to do the operations with a precision of 2-( m+ °^ lo s J')) in order for the final 
coefficients to be with precision 2~ m . 

Note that in the case when / is a finite degree polynomial, we can evaluate the coefficients 
of its powers using the multinomial formula, and with no need for the numerical iterative 
computation described above. 

Once we have the coefficient cp of in f{z) 1 , we are able to write 

=a 3 + at l) cf + at^cf + ... + df^^. (7) 

We already have all the parameters in (JJJ), as numbers or explicit polynomials in n of degree 
at most j — 2, except for . It is easy to see that Cj = 1. Thus, we obtain an explicit 
recurrence, that connects aj 1 ^ with CLf _1 \ and yields 



af = aj + J2 («i + afcf + afcf ] + ... + af^c^ 

i=l 

Thus, is given by a polynomial of degree at most j — 1 in n. The coefficients can 
be computed very efficiently (see [GKP94J for more information on how to compute the 
sum J2i=i^)- The precision bit loss in this process is also limited to 0(k\og 2 k) bits. The 
coefficients of the polynomial are precisely the information we are looking for to complete 
the proof of the lemma. I 

6 Computing a "long" iteration 

We can now apply the results of sections |H and El to prove the following lemma. 
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Lemma 8 Suppose f(z) = z + a r z r+1 + a r+ iz r+2 + . . . is given by its power series with some 
positive radius of convergence R. Then there is an easily computable number C such that 
if \z\ < i < R, we can compute the £ = ^-\-th iterate of z and its derivative ^jj(z) with 
precision 2~ s in time polynomial in s and logm. 

The loss of precision from z to f e (z) can be bounded to a constant number of bits. The 
loss of precision from z to ^-(z) is 0(—\og\z\) bits. 

Proof: We begin similarly to the discussion in the beginning of section As before, 
it is easy to compute a power of two, A = 2 a such that |a$| < A 1 " 1 for all i. Again, let 
g{z) = Af(z/A). Then all the coefficients of g are bounded by 1 in absolute value. Write 

g i (z) = z + b®z r+1 + b ( $ 1 z r+2 + ... 

g is dominated by the series z + z r+1 + z r+2 + . . .. Thus, using lemma El we conclude that 
| bf>\ < {atfl r for some simple, computable a. If we write 

f(z) = z + a^z r+1 + + . . . = jg £ (Az), (8) 

we see that = A k o£\ and \a^ \ < (aA r £) k / r . Considering that £ = and \z\ < ^, we 
obtain 



k 

\Cm r+1 J " V C 1 ^') ' 

Choose C > 2 r aA r . Then |a^z fc+1 | < 2~ fc , and it suffices to consider the first s + 2 terms 
of the series (jHJ) to obtain the desired iteration with a 2~ s precision (all later terms become 
negligible). 

All we have to do now is to compute af> , <4+i5 • • •> a i+2 with precision 2~( s+e M). We do 
it by computing their coefficients from (jHJ) with precision 2~^ + °( slog£ ^, which can be done 
in time polynomial in s and \og£ = B(logm) by lemma 

To compute the derivative of the £-th iteration, write 

^-( z ) = 1 + ( r + l) a <V + (r + 2)a%z r+l + ... 

then 

/ n A r m r \ k/ r 

\{k + l)afz k \ <(k + l) {^r) < (k + 1) • 2- fe < 2^, 

for sufficiently large k. Hence it suffices to consider the first 2s + 2 terms of the series (jHJ) 
to obtain the desired iteration with a 2~ s precision (all later terms become negligible). We 
compute a£+i, • ■ ., with precision 2~( s+e ^\ which again can be done in time 

polynomial in s and \og£ = G(logm) by lemma 

The loss of precision can be kept to a constant number of bits for z by the constant 
bound we have on the first derivative of f (z) around z. The second derivative is bounded 
by O (l/\z\) around z, and the precision loss can be kept to 0(— log \z\) bits, which is fine, 
as long as \z\ is not too small. H 
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7 Computing parabolic Julia sets in polynomial time 



In this section we put the pieces together to give a poly-time algorithm for computing 
parabolic Julia sets. For the rest of the section fix r(x) to be a rational function on <D, and 
denote its Julia set by J r . We consider J r first as a subset of the Riemann sphere C. Using 
a stereographic projection n : C — {00} — > C, for any compact CcC such that 00 C, it 
is easy to see that computing J r on C is exactly as hard as computing tt( J r ) in C. Hence, if 
00 is not in J r , it suffices to compute it in some bounded region B(0, R) in C. 

If 00 e J r , then it is obviously impossible to "draw" it on the plane. We can still "draw" 
it on the Riemann sphere, and hence on any bounded region of <D. We take a Mobius 
transformation T such that for the conjugation r' = T o r o T -1 , J r .i is obtained from J r by 
a rotation of the Riemann sphere. In this way, drawing J r i on C is as easy (or as difficult) 
as drawing J r . We can choose T so that 00 ^ J r i, and then it suffices to draw J r > on some 
bounded region of C 

From now on, we assume that 00 ^ J r , and that we have some constant B such that 
J r C B(0, B) C C. We are trying to compute J r on this bounded region. 

7.1 Preliminaries — the nonuniform information we will need 

Below we list the information the algorithm will use to compute J r efficiently with an arbi- 
trarily high precision. We will need the following ingredients: 

1. A list of periods v 1, v 2, . . . , v k for all the parabolic orbits. Consider the iteration r v (z) 
of r(z) for v = LCM(t>i,f 2 , . . . ,Vk)- It has only simple parabolic points (no orbits). 
From now on, we replace r(z) with r v (z). We can do it, because J r v = J r for all v. 
We can multiply v by some other factor, so that the derivative dr ^ at each parabolic 
point is 1, and not some other root of unity. 

2. Information that would allow us to identify the parabolic points, and information about 
them. For each parabolic point p we would like to know an approximation q for p that 
would allow us to compute p in poly-time using Newton's method. Near p, r(z) can 
be written as 

r(z) =p+(z-p) + a u (z - p) u+l + a u+1 (z - p) u+2 + ... 
for some integer u. We would like to know this number. 

3. An open set U such that 

• U DV := r~\U), 

• all the parabolic points are in dU, 

• J r C U, 

• all the critical (and hence also all the postcritical) points of r(z) lie outside U, 
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• all the poles and their neighborhoods lie outside U, 

• moreover, if we denote the postcritical set by P r , then for any u G r(U), d(u, P r ) > 
32 ■ d(u, J r ), and 

• outside any e-neighborhood of the parabolic points and a finite number of their 
pre-images, the distance between dll and dV is bounded from below by some 
positive 5. 

U is given in the form U = r~ u (U), where U is some explicit semi- algebraic set. Thus 
queries about membership in U and V can be computed efficiently with an arbitrarily 
high precision, at least outside some small region around the parabolic points and their 
preimages up to order u. We will show that such a U exists, and how to compute it 
from some basic combinatorial information in section 17731 

4. For each parabolic point p, there is a small neighborhood E p of p in which lemma |H1 
applies for computing a long iteration of r. We would like to have two sets E\ and E 2 
around the parabolic points and their pre-images of order up to u such that 

• E 1 C E 2 , 

• for a given point z, it takes constant time to decide whether z G E\, or z £ E 2 , 

• for each z G E 2 , there is a parabolic point p such that w = r u (z) G E p , 

• OUndV C E u and 

• we have a positive d such that for any two points X\ G V, x 2 U outside of E\, 
\xi — x 2 \ > d. 

5. The set dll D dV consists of pre-parabolic points q, i.e. points such that r u (q) is 
parabolic for some fixed q. The repelling directions and their pre-images belong to V. 
There is an angle a such that all the points in E 2 that form an angle < a with one 
of the repelling directions, or their preimages belong to V. If necessary, we can make 
E 2 smaller. We denote the subset in E 2 of points that make an angle of < a/2 with a 
repelling direction or its pre-image by A±, and the points that make an angle of < a 
with a repelling direction or its pre-image by A 2 . We can choose a as small as we want. 
We have the following properties. 

• A x c A 2 c V, 

• if z is given within an error of < \z — q\ 2 , near a pre-parabolic point q G r~ u (p), 
we can tell if z G A\ or z A 2 , 

• for any p, and for any w G A 2 D E p , \r'(w)\ > 1. This is true for a sufficiently 
small a. 

6. Consider the Poincare metric defined on the hyperbolic set U. Denote its density by 
djj. We have the following theorem, known as Pick's theorem (see [MilOOj for a proof). 



13 



Theorem 9 (Theorem of Pick) Let S and T be two hyperbolic subsets of C If 
f : S — > T is a holomorphic map, then exactly one of the following three statements is 
valid: 

(a) f is a conformal isomorphism from S onto T, and maps S with its Poincare 
metric isometrically to T with its Poincare metric. 

(b) f is a covering map but is not one-to-one. In this case, it is locally but not globally 
a Poincare isometry. Every smooth path P : [0, 1] — > S of arclength I in S maps 
to a smooth path f o P of the same length I in T . 

(c) In all other cases, f is a strict contraction with respect to the Poincare metrics 
on the image and preimage. 

Let dy be the density of the Poincare metric defined on V. By the construction, V 
contains no critical points, and so r : V — > U is a covering map and by theorem |U] it is 
a local isometry. That is, for any z G V, 

d v {z)=d u {r{z))-\r'{z)\. (9) 

On the other hand, the embedding i : V <— »• U is not a covering map, hence it is strictly 
contracting in the Poincare metric. Thus for any z G V we have dy(z) > djj{z). 
Together with Q, this implies 

d u {r{z))-\r'{z)\ = d v {z)>d u {z) (10) 

for all z & V. In particular, if we consider only z G V — Ei, then z is in some compact 
domain bounded away from the boundary of U, hence the ratio dy(z)/djj(z) is always 
positive (maybe oo), and it has a minimum c > 1. We would like to have this c as 
part of the nonuniform information. With this c we have dy{z) > c ■ du{z) for all 
z G V — Ei, and (fTUj) becomes 

du(r(z)) ■ \r\z)\ = d v (z) > c ■ d v {z) (11) 

for all z G V — Ei. Moreover, we can choose a slightly smaller c such that (fTTj) holds 
for any point z on any path p from w G V — E 1 to J r such that L d[J {p) < 2dd u (w, J r ). 
This is true since the lengths of such paths can be uniformly bounded, and thus it 
cannot get too close to the points of dU PI dV . We would like to have the value of c 
(or some rational estimate Co, 1 < Co < c). 

7. Since the postcritical points are outside U, and the parabolic points cannot be critical, 
we can have a constant d > such that |r'(;z)| > d for all z G r(U) U E2 (if necessary, 
we can choose a smaller E2. 

8. Finally, we need an efficient procedure to estimate the distance from J r for all points 
that are not too close to it. More specifically, for any point z outside of VUEi, there is 
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an "estimator" that provides the distance d(z, J r ) within a multiplicative error factor 
of 2. This can be done since J r C V and the distance d(d(V U Ex), J r ) is bounded 
from below by a constant. Hence a fixed-precision image of J r suffices to make such 
an estimation. 

The situation is somewhat different if z G E 2 . In this case, we know that either z or 
r k (z) for some bounded k is close to a parabolic point p. We know that in some small 
neighborhood of p, J r looks like I lines at angle 27r/£ from each other leaving p (see 
[MilOOj, Chapter 10 for more details). Denote this set by L p . In general, if z is very 
close to the Julia set, d(z,L p ) can be (multiplicatively) very different from d(z,J r ). 
However, if we stay away from the repelling directions (and hence from L p and J r ), 
these two quantities are actually similar. More precisely, for any a (and in particular 
for a mentioned in the definition of Ax), in a small neighborhood of p, d(z, L p ) is within 
a factor of 2 from d(z, J r ) for every z that has an angle of at least a/2 with each of the 
repelling directions. The same holds for the first u pre-images of the parabolic points. 
We can take E 2 (and hence Ax) to be sufficiently small such that this property holds 
within E 2 — Ax. 

7.2 Algorithm outline and analysis 

The goal of the algorithm is to compute a function from the family 



To do this, we estimate d(z, J r ) up to a multiplicative constant, assuming that d(z, J r ) > 2~ n . 
If the assumption does not hold, the algorithm always outputs 1 (or successfully estimates 
the distance). 

The algorithm outline is as follows. 

1. w <— z; steps <— 0; 

2. ifw^VUEx, output 0; 

3. estimate the maximum number iV = 0(n) of steps outside E\ we would need; 

4. iterate the point w as follows: 

5. set derivative counter D <— 1; the cumulative derivative estimation should be bounded 
between the derivative and twice the derivative at all steps; 

6. if w i VVJEx: 

(a) estimate e, a 2-approximation of d(w, J r ): e < d(w, J r ) < 2e; 

(b) output 1 if % < 8 • 2' n , and if % > 16 • 2~ n ; 




1, if d(z, J r ) < 2' n 

0, if d{z, jj) > 256 • 2- 
or 1, otherwise 



n 
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7. ifweU-Er. 

(a) D^D - \r'(w)\ 

(b) w <— r(w); 

(c) steps <— steps + 1; 

(d) if steps > N, output 1; 

8. if w G E 2 in the region of some preparabolic q G r~ u (p), and |w — g| < 2~@ n {(3 - a 
constant to be determined): 

(a) output 1; 

9. if w G E 2 PI v4 2 , and it is not in the neighborhood of any parabolic point: 

(a) D^D- \r'(w)\ 

(b) u> <— r(u>); 

10. if iu G £ 2 - Ai: 

(a) estimate e, a 2-approximation of d(w, J r ): e < d(w, J r ) < 2e; 

(b) output 1 if % < 8 • 2" n , and if § > 16 • 2~ n ; 

11. if w G -Ep fl A 2 for some p: 

(a) make a long iteration ?/ = r v (w); 

(b) if ?/ G -Ep, but escapes A 2 : 

i. do binary search to find the smallest u < v such that r u (w) is in A 2 — A±; 

ii. w <— r u (w) 

iii. go to step [TU| 

(c) else, w <— r^(w); 

(d) D <- D ■ 



dr v (w) 
dw 



The algorithm performs the operations with 0(n 2 ) bits of precision. First note that steps 
imrm cover all the possibilities for w. If two or more of the possibilities intersect, it does not 
matter which one to choose. 

We first show that 

Claim 10 Step&in the algorithm is possible. 

Proof: Let z be some point in V outside Ex. Let p be the shortest path in the Poincare 
metric djj from p(0) = f(z) to p(l) G J r . r : V — > U is a covering map, and p can be raised 
to a path p in V such that p(0) = z and p(l) G J r (by the invariance of J r under r). There 
are two possibilities: 
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1. L du (p) > 2d du (z, J r ). We know that r expands the Poincare metric djj, and hence 



d du {r(z), J r ) = L du (p) > L du (p) > 2d du (z, J r ) > c ■ d du (z, J r ). 



2. Otherwise, by property H of c from section I7~T1 r is expanding by a factor of c along 
the entire path p. Hence 



This shows that every step [7| multiplies the Poincare distance between w and J r by a factor 
of at least c. Other steps do not decrease it. If initially the Euclidean distance d(w, J r ) is at 
least 2 _n , then the Poincare distance is at least C ■ 2~ n for some constant C. For every point 
in V — Ei this distance is bounded from above, hence it will take at most log c C _1 -2 n = 0(n) 



Claim 11 At any stage of the algorithm, w G r(U). 

Proof: If on the first iteration we do not exit on step |H1 or EH then we must have 
w G A 2 UV = V C U before the first iteration. From here, we prove the claim by induction. 

Suppose the algorithm is running after i iterations. Denote the current value of w by w' 
and the value after the iteration by w" (w" = w' if the algorithm terminates). We assume 
that w' G r(U). If the algorithm executes steps [7| or El then by the conditions w' G U, and 
w" = r(w') G r{U). Quitting on steps IH1 and ITU1 does not affect the value of w. Step [TT] only 
runs to keep w" in A 2 C V. Thus w" G r{U) in this case as well. I 

The following is a classical theorem in complex analysis. 

Theorem 12 Koebe's 1/4 Theorem Suppose (f) : (Si,Si) — ► (5*2,52) is a holomorphic 
bijection between two simply connected subsets of C, Si and S 2 - Let r\ be the inner radius 
of Si around Si, and r 2 be the inner radius of S 2 around s 2 . Then the following inequality 
holds: 



We apply Koebe's theorem to prove the following lemma. 

Lemma 13 Suppose that z G r(U) and the algorithm terminates with w = r k (z), D < 
dr d ^ < 2D, and e < d(w, J r ) < 2e for some distance estimate e. Then 



d du (r(z), J r ) = L du {p) >c-L du {p)>c - d du (z, J r ). 



(12) 



steps [7| for the orbit of z to escape. 



r 2 > jri ■ \(j>'(si)\. 



< d(z, Jr) < 

8D ~ v ' ' ~ 



16e 
~D 



(13) 
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Proof: Denote a = dl ' ^ and s = d(w, J r ). Then since w E r{U) the distance d(w,P r ) 
from the postcritical set is at least 32s. We consider orbit z, r(z), r 2 (z), . . ., r k (z) = w. 
Let St = B(w,32s). Consider the preimage Sk-i of Sk under the branch of r that takes 
r k ~ 1 (z) to r k (z). It is uniquely defined since Sk contains no postcritical points. The mapping 
r : Sk-i — > Sk is a one-to-one conformal mapping. We can continue this process to obtain 
a one-to-one conformal branch r h : (5*0,-2) — * (Sk,w). By Koebe's theorem the image of 
B(w,s) under the inverse of this mapping must contain a ball of radius at least around 
z. By the invariance of J r , this ball contains no points from J r . Hence 

7/ T \ S 6 6 

a(z, J r ) > — > = • 

1 ' r > - 4a - 4 • 2D 8D 

Also by Koebe's theorem, So contains the ball B\ of radius = — around z. The image 
r k {B\) must contain a ball of radius — - 1 = 2s around w. Hence r k (Bi), and also B\ contain 
points from J r . So 

d(z,J r )<r(B x ) = — < — = —. 



Claim 14 If the algorithm terminates at step\U or MUL it outputs a valid 



answer. 



Proof: The variables e and D in these cases satisfy the conditions of lemma El If 
d(z, J r ) < 2~ n , then < d(z, J r ) < 2~ n , and ^ < 8 ■ 2 _n , so the algorithm outputs 1. If 
d(z, J r ) > 256 ■ 2~ n , then ^ > d(z, J r ) > 256 ■ 2~ n . Hence % > 16 ■ 2" n , and the algorithm 
outputs 0. I 



Claim 15 SYej>^| is executed at most N + 1 times, and if it outputs 1, it is a valid answer. 

Proof: This follows from the definition of the number N, the existence and computability 
of which has been established in claim M 

Claim 16 Step\^is executed at most a constant number of times between two executions of 
step^ 

Proof: This is true because every point in E2 is either in the neighborhood of a parabolic 
point, or a pre-image of order at most u of a parabolic point. Hence we can have at most u 
iterations of step M before an iteration with step [7| or [TT] being executed. A series of step ITT1 
iterations ends with a termination or with a step [7| iteration before another step 01 iteration. 

■ 

Claim 17 After j iterations, out of which i are of step [7| or \Q \D\ > d l . In particular, by 
claimsEB andUM we always have \D\ > d°^ = d°^ = 2~ 0(n \ 

Proof: By property each step [3 and M contribute at most d to D. Steps IB| IH1 and fTUl 
terminate the algorithm. In step ^3 we do the long iteration only until the angle between 
a repelling direction and w exceeds a. Until that moment, by property El in section 17.11 
|r'(w)| > 1, hence step [TT] does not decrease D. M 



18 



Claim 18 Step\& outputs a valid answer for some constant (3, for sufficiently large n. 



Proof: This is clearly true if w = z it the 0-th iteration. Otherwise, it is obvious that the 
previous step could not have been a step[TTJ hence it must have been either a step [7| or El In 
either case, during the previous iteration the value of w was some w' such that r(w') = w. 
By claim ITT1 w' G r(U). 

Denote the preparabolic point near w by q. Then w' is a preparabolic point q' of order 
at most u + 1. Since the parabolic points are not postcritical, there is some 7 such that a 
27 neighborhood of any preparabolic point of order < u + 1 is mapped by r in a one-to-one 
fashion with no critical points. > d, since q' G J r C r(U), and by Koebe's theorem, 

r(i?(g',7)) contains a ball of radius d ■ 7/4 around q. We can take (3 sufficiently large, so 
that 2~ Pn < d ■ 7/4 for all n. 

Suppose the algorithm exits on step [HI Denote d! = d(w, J r ) < \w — q\ < 2~ l3n < d ■ 7/4. 
Consider the one-to-one restriction f of r to B{q',2 r j). w is in the d ■ 7/4-neighborhood 
of q, and hence w' is in the 7-neighborhood of q'. Denote e = d(w',J r ). e < 7, since 
q' G J r - Consider the set r(B(w', e/2)). It contains no points of J r , and by Koebe's theorem 
it contains the ball B(w,e\r'(w')\/8). Hence d(w, J T ) > e\r'(w')\/8 > e ■ d/8. On the other 
hand d(w,J r ) < 2~ /3n . Combining these inequalities we obtain 2~ /3n > e ■ d/8, and hence 
e < 8 ■ 2~ l3n /d. 

By claim El we always have D > 2~ T ' n for some constant 77. w' G r(U), and by lemma 
El with w' we have 

d( z , J r )<^< ^'^^ = 2 7 ~ f5n ~ lo ^ n < 2-", 

if we take j3 > 8 — \ogd + 77. This is the value of (3 we should take. H 

Claim 19 Step[Tl\is executed at mostO{n) number of times between two executions of step 


Proof: According to lemma |H1 (the long iteration lemma), we can make a step Q(m r ) 
iterations if \z\ < — . For a sufficiently small a, and close enough to the parabolic point, if 
\z\ > 2"-:, we have 



\ r m r /C 



(z)\ = \z\+n(m r \z\ r+1 /C) = \z\+Q ( ^ ) = \ Z \ + Q ( } ) > \z\-(l + 5), 



for some 5 > 0. If the algorithm did not terminate at step El \z\ > 2 _/3n , and it will take 
0(n) long iterations to escape the neighborhood of the parabolic point and either terminate 
or reach a step I 

It follows from the claims that 

1. The algorithm terminates after 0(n 2 ) iterations. 

2. When it terminates, it outputs a valid answer. 

This shows that the algorithm is polynomial and correct. 
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7.3 Uniformizing the construction 

In this section we show how to uniformize the construction. In other words, we are trying 
to construct one machine computing J r for the biggest possible family of parabolic r's. As 
has been mentioned in section 12.21 the output of the machine varies continuously in the 
Hausdorff metric with the input coefficients. The map J : c i— > J z 2 +C is discontinuous at the 
parabolic point c = 1/4 (see |Dou94j ) . Thus, we cannot expect one machine to compute all 
hyperbolic and parabolic sets even in the quadratic case. 

Despite the big number of different parameters that were mentioned as pre-requisites in 
section I7TT1 we will argue that all the information can be derived from some basic information 
about the number and periods of the parabolic points. 

First, we prove the following. 

Claim 20 Given the information on the parabolic orbits, we can extract the information on 
the attracting orbits ourselves. 

Proof: The immediate basin of each attracting periodic orbit contains at least one critical 
point (e.g. Theorem 8.6 in |MilOO| ). On the other hand, in our case every critical point 
converges either to an attracting or to a parabolic orbit. We proceed as follows. Iterate 
each critical point until we know to which orbit it converges. If it converges to an attracting 
orbit, we will eventually know it. Continue this process until the convergence of all critical 
points is accounted for. H 
Probably the most interesting part is computing the set U. So far we haven't even shown 
that such a U exists. To compute U, we start with a set U defined as follows. Around each 
attracting orbit we take a small ball in the basin of attraction. Denote the union of these 
balls by A. 

Around each parabolic point p, for any attracting direction d, consider a small "diamond" 
shaped region P around d such that: 

• r(P) C P with P n r(P) = {p}, and 

• the angle of P at p is at least || of the angle between the two repelling directions. 

The edges of P are chosen so that points would not escape it under r. See figure El for an 
illustration. This is possible by the basic properties of the series expansion of r near p. 
Denote the union of these "diamons" P by P. Define 

u = c -{AuP). 

We know that there is an iteration q such that 

1. for all critical points c, r 9_1 (c) is outside U, and 

2. r 9_1 (oo) is outside U. 
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Figure 3: The region P, here there are 3 repelling directions 

This is true, since the orbits of all the critical points eventually converge either to an at- 
tracting or a parabolic orbit, and by our assumption so does the orbit of oo. We have the 
following claims. 

Claim 21 V = r~ l (U) C U , and dU n dV = {the set of parabolic points}. 

Proof: This follows immediately from the definition of U . I 

Claim 22 

oo 

(1 r- n (U) = J r . (14) 

n=0 

Proof: The orbit of any z G J r always stays in J r C U, hence such a z is in the intersection 
above. 

The orbit of any w ^ J r eventually converges to either an attracting or a parabolic orbit, 
and thus escapes U. This means that r k (w) ^ U for some k, and w ^ r~ k (U). So w is not 
in the intersection in this case. H 

Claim 23 For sufficiently large q, U = r~ q (U) satisfies the conditions of part&in section 

123 
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Proof: The first five conditions are satisfied automatically by the definition of U. The 
last condition follows from the fact that for any q, d(r~ q (U)) D d(r~ Q (V)) consists of the 
parabolic points and their pre-images up to order q. 

The hardest condition to satisfy is the sixth one. Namely, we want to have for any 
u G r(U), d(u, P r ) > 32 • d(u, J r ). Here P r denotes the postcritical set of r. 

Let N e be an e-neighborhood of the parabolic points, and M £ - an e-neighborhood of the 
attracting orbit points. We know that for any e only finitely many points of P r lie outside 
N £ UM £ . For a sufficiently small e all the points in P r C\N e lie in a small angular neighborhood 
of the attracting directions. 

Denote by 7 the angle between two adjacent repelling directions. By the definition of 
P, for any q, all the points in r~ q (U) are in a ^|g7-neigborhood of the attracting direction. 
Thus the condition is satisfied in N £ . 

Outside of N £ U M E , there are only finitely many points of P r , hence there is a minimum 
d of their distances from J r . This minimum also exists for points in M £ , since attracting 
orbits are bounded away from J r . By (fT4"j) and compactness, for a sufficiently large q, r~ q (U) 
is in the ^-neighborhood of J r , and the condition is satisfied outside of N £ . M 

Claim 24 For any constant c, we can produce a 2~ c -precise image of J r . 

Proof: This can be done using a procedure described in BBY04J, Theorem 1.2. Note 
that since c does not depend on n, the running time of this procedure will not depend on n 
as well. M 
In particular, claim |^] immediately allows us to establish property |H1 in section 17.11 

Claim 25 The q (and hence U) from claim can be computed from the basic information 
about the parabolic points. 

Proof: The proof of claim|2Slis constructive, except for the argument outside of N e , which 
uses compactness. There are only finitely many points of P r outside of N £ U M £ , which can be 
easily computed. The distance from J r to the points of M £ can also be easily bounded from 
below. We can find the desired value of q by computing a sufficiently good approximation 
of J r , which is done by claim The precision with which we will have to perform this 
computation depends on r but not on n. I 

The other parts of the construction are easily seen to be uniformizable. 
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